Aim and Study Area

From a 1m DSM in Tiraumea, New Zealand, terrain derivatives will be computed with the following code. The aim is to detect Earthflows, possibly through machine learning, where a model would be fed with the terrain derivatives, and optical information.

Since the data is large, two initial study areas were proposed for analysis:

# Call mapping library and set options
library(mapview)
mapviewOptions(basemaps = 'OpenStreetMap', console = F, verbose = F, leafletWidth = 800)
library(sf, quietly = T, verbose = F, warn.conflicts = F)
tiraumea = st_read('study_area/Tiraumea.shp', quiet = T) 
sa = st_read('study_area/study_area_outlines.shp', quiet = T)
mapview(list(tiraumea,sa))

Workflow

Initially, we will look at Study Area 1. I will load the DEM data into R with the help of the raster package. The DEM sinks were previously filled with this example code:

rsaga.fill.sinks(in.dem = "dem.sgrd", out.dem = "dem_filled", method = "wang.liu.2006", env = env)

Then it was cropped to the study area (in QGIS).

Now for the resulting DEM:

# I use 'raster` at this point, since for some reason writing a `stars` object into a SAGA compatible format is not working. 
library(raster, quietly = T)
dem = raster('input_dem/dsm_filled_sa1.tif') 
crs(dem) = '+init=epsg:2193'
mapview(dem)

I save it as a SAGA compatible format to be able to call it into the RSAGA functions.

writeRaster(dem, 'input_dem/dsm_sa1.sgrd', format = 'SAGA', overwrite = T, NAflag = 0, prj = T)

Generating Terrain Derivatives

First I set up the RSAGA environment:

library(RSAGA)
#library(link2GI)
#saga = linkSAGA()
#saga # Copy the results of this to the env setting
env = rsaga.env(
  path = "C:\\Users\\b1066081\\Desktop\\saga-7.5.0_x64",
  modules = "C:\\Users\\b1066081\\Desktop\\saga-7.5.0_x64\\tools"
  # cmd =  "C:\\Users\\b1066081\\Desktop\\saga-7.5.0_x64\\saga_cmd.exe"
)
Verify specified path to SAGA command line program...
Verify specified path to SAGA modules...
Done

And now, compute some derivatives for the different modules.

Morphometry module

Slope, Aspect, Curvature module

With this I will compute overall 8 outputs:

  • Slope (slope) in degrees
  • General Curvature (cgene)
  • Plan Curvature (cplan)
  • Profile Curvature (cprof)

I will also compute variables to keep with the workflow of Eisank et al. (2014), but I will not add the multiscale parameter, since I do not find how this was done for that.

  • Crossectional curvature (ccros)
  • Longitudinal curvature (clong)
  • Max. curvature (cmaxi)
  • Min. curvature (cmini)
rsaga.slope.asp.curv(
  in.dem = "input_dem/dsm_sa1.sgrd", 
  out.slope = "out_products/slope_sa1.sgrd",
  out.cgene = "out_products/cgene_sa1.sgrd",
  out.cplan = "out_products/cplan_sa1.sgrd",
  out.cprof = "out_products/cprof_sa1.sgrd",
  out.ccros = "out_products/ccros_sa1.sgrd",
  out.clong = "out_products/clong_sa1.sgrd",
  out.cmaxi = "out_products/cmaxi_sa1.sgrd",
  out.cmini = "out_products/cmini_sa1.sgrd", 
  unit.slope = "degrees", env = env,
  method = "poly2zevenbergen"
)

Relative heights and slope positions

From this modules on, I will need to find out the usage first, because there are no shortcuts for them. I will always refer to this table for the morphometry module to get the respective code:

rsaga.get.modules('ta_morphometry', env = env) %>% knitr::kable()
code name interactive
0 Slope, Aspect, Curvature FALSE
1 Convergence Index FALSE
2 Convergence Index (Search Radius) FALSE
3 Surface Specific Points FALSE
4 Curvature Classification FALSE
5 Hypsometry FALSE
6 Real Surface Area FALSE
7 Morphometric Protection Index FALSE
8 Multiresolution Index of Valley Bottom Flatness (MRVBF) FALSE
9 Downslope Distance Gradient FALSE
10 Mass Balance Index FALSE
11 Effective Air Flow Heights FALSE
12 Diurnal Anisotropic Heat FALSE
13 Land Surface Temperature FALSE
14 Relative Heights and Slope Positions FALSE
15 Wind Effect (Windward / Leeward Index) FALSE
16 Terrain Ruggedness Index (TRI) FALSE
17 Vector Ruggedness Measure (VRM) FALSE
18 Topographic Position Index (TPI) FALSE
19 TPI Based Landform Classification FALSE
20 Terrain Surface Texture FALSE
21 Terrain Surface Convexity FALSE
22 Terrain Surface Classification (Iwahashi and Pike) FALSE
23 Morphometric Features FALSE
24 Valley and Ridge Detection (Top Hat Approach) FALSE
25 Fuzzy Landform Element Classification FALSE
26 Upslope and Downslope Curvature FALSE
27 Wind Exposition Index FALSE
28 Multi-Scale Topographic Position Index (TPI) FALSE
29 Wind Shelter Index FALSE 
rsaga.get.usage('ta_morphometry', 14, env = env)

From here I will get 5 outputs, and will work with default W, T, E.

  • Slope height (slhgt)
  • Valley depth (vldpt)
  • Normalized height (nrhgt)
  • Standardized height (sthgt)
  • Mid-slope position (mdslp)
rsaga.geoprocessor(
  'ta_morphometry', 14, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    HO='out_products/slhgt_sa1.sgrd', 
    HU='out_products/vldpt_sa1.sgrd',
    NH='out_products/nrhgt_sa1.sgrd',
    SH='out_products/sthgt_sa1.sgrd',
    MS='out_products/mdslp_sa1.sgrd'
  ),
  env = env
)

Downslope distance gradient

This tool has its own module, and will return one outcome, coded ‘ddgrd’. We get the usage:

rsaga.get.usage('ta_morphometry', 9, env = env)

And call it, with default distance (10) and gradient in degrees:

rsaga.geoprocessor(
  'ta_morphometry', 9, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    GRADIENT='out_products/ddgrd_sa1.sgrd', 
    OUTPUT = 2
  ),
  env = env
)

Mass Balance Index

This tool has its own module, and will return one outcome, coded ‘mbidx’. We get the usage:

rsaga.get.usage('ta_morphometry', 10, env = env)

And compute with all the defaults:

rsaga.geoprocessor(
  'ta_morphometry', 10, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    MBI='out_products/mbidx_sa1.sgrd' 
  ),
  env = env
)

Topographic Position Index

This tool has its own module, and will return one outcome, coded ‘tpidx’. We get the usage:

rsaga.get.usage('ta_morphometry', 18, env = env)

And compute with all the defaults, min and max are definitely needed and vary depending on the raster resolution. If not used properly computation will be exhaustive.

rsaga.geoprocessor(
  'ta_morphometry', 18, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    TPI='out_products/tpidx_sa1.sgrd',
    RADIUS_MIN=0,
    RADIUS_MAX=20
  ),
  env = env
)

Convergence Index

This tool has its own module, and will return one outcome, coded ‘cvidx’. We get the usage:

rsaga.get.usage('ta_morphometry', 1, env = env)

And compute with method gradient and neighbours 3x3:

rsaga.geoprocessor(
  'ta_morphometry', 1, 
  list(
    ELEVATION='input_dem/dsm_sa1.sgrd', 
    RESULT='out_products/cvidx_sa1.sgrd', 
    METHOD=1,
    NEIGHBOURS=1
  ),
  env = env
)

Terrain Roughness Index

This tool has its own module, and will return one outcome, coded ‘tridx’. We get the usage:

rsaga.get.usage('ta_morphometry', 16, env = env)

And compute with all the defaults:

rsaga.geoprocessor(
  'ta_morphometry', 16, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    TRI='out_products/tridx_sa1.sgrd'
  ),
  env = env
)

Texture

This tool has its own module, and will return one outcome, coded ‘textu’. We get the usage:

rsaga.get.usage('ta_morphometry', 20, env = env)

And compute with all the defaults:

rsaga.geoprocessor(
  'ta_morphometry', 20, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    TEXTURE='out_products/textu_sa1.sgrd'
  ),
  env = env
)

Hydrology module

Again I will access the codes for the tools to get the usage:

rsaga.get.modules('ta_hydrology', env = env) %>% knitr::kable()
code name interactive
0 Flow Accumulation (Top-Down) FALSE
1 Flow Accumulation (Recursive) FALSE
2 Flow Accumulation (Flow Tracing) FALSE
4 Upslope Area FALSE
6 Flow Path Length FALSE
7 Slope Length FALSE
10 Cell Balance FALSE
13 Edge Contamination FALSE
15 SAGA Wetness Index FALSE
16 Lake Flood FALSE
18 Flow Accumulation (Mass-Flux Method) FALSE
19 Flow Width and Specific Catchment Area FALSE
20 Topographic Wetness Index (TWI) FALSE
21 Stream Power Index FALSE
22 LS Factor FALSE
23 Melton Ruggedness Number FALSE
24 TCI Low FALSE
25 LS-Factor, Field Based FALSE
26 Slope Limited Flow Accumulation FALSE
27 Maximum Flow Path Length FALSE
28 Flow between fields FALSE
29 Flow Accumulation (Parallelizable) FALSE

Melton Ruggedness Index

This tool has its own module, and will return one outcome, coded ‘mridx’. We get the usage:

rsaga.get.usage('ta_hydrology', 23, env = env)

And compute with all the defaults:

rsaga.geoprocessor(
  'ta_hydrology', 23, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    MRN='out_products/mridx_sa1.sgrd'
  ),
  env = env
)

LS Factor

This particular tool I had to compute using the GUI for SAGA, because RSAGA does not let me use the LS Factor (one step) library that I need to do it right. In the GUI I only need to give the DEM. It is saved as ‘lsfct’.

SAGA Wetness Index

This tool has its own module, and will return one outcome, coded ‘sagaw’. We get the usage:

rsaga.get.usage('ta_hydrology', 15, env = env)

And compute with defaults:

rsaga.geoprocessor(
  'ta_hydrology', 15, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    TWI='out_products/sagaw_sa1.sgrd'
  ),
  env = env
)

Slope Length

This tool has its own module, and will return one outcome, coded ‘sllgt’. We get the usage:

rsaga.get.usage('ta_hydrology', 7, env = env)

And compute with defaults:

rsaga.geoprocessor(
  'ta_hydrology', 7, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    LENGTH='out_products/sllgt_sa1.sgrd'
  ),
  env = env
)

Stream Power Index

To compute this index, some previous steps are required. The version incompatibility of SAGA 7.5 and RSAGA does not support this computation(?), at least in an intuitive way. So I compute on the SAGA GUI.

The steps I took were:

  1. Calculate Flow Accumulation (One Step) on the ta_hidrology module. The only input is the DEM. The intermediate output needed is the Specific Catchment Area. I saved this on the int_products directory as ‘spcar_sa1’. Specific parameters are:
    • Preprocessing: Fill Sinks (Wang & Liu)
    • Flow Routing: Deterministic 8
  2. Compute Stream Power Index on the ta_hidrology module. Inputs are the Slope and Specific Catchment Area. The result is saved to out_products as ‘spidx’.

Terrain Wetness Index

This particular tool I had to compute using the GUI for SAGA, because RSAGA does not let me use the TWI (one step) library that I need to do it right. In the GUI I only need to give the DEM. I selected Deterministic 8 for the Flow Distribution. It is coded ‘twidx’.

Lightning module

Again I will access the codes for the tools to get the usage:

rsaga.get.modules('ta_lighting', env = env) %>% knitr::kable()
code name interactive
0 Analytical Hillshading FALSE
2 Potential Incoming Solar Radiation FALSE
3 Sky View Factor FALSE
4 Topographic Correction FALSE
5 Topographic Openness FALSE
6 Visibility (points) FALSE
7 Potential Annual Insolation FALSE
8 Geomorphons FALSE

Sky View Factor **

For this and the next tool, a multiscale parameter is chosen to speed up the computation. From this tool we can get two derivatives:

  • Sky view factor (svfct) - [not simplified]
  • Visible sky (visky)
rsaga.get.usage('ta_lighting', 3, env = env)

And compute with defaults:

rsaga.geoprocessor(
  'ta_lighting', 3, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    VISIBLE='out_products/visky_sa1.sgrd',
    SVF='out_products/svfct_sa1.sgrd',
    METHOD=1,
    DLEVEL=3
  ),
  env = env
)

Topographic Openness **

From this tool we can get two derivatives:

  • Positive openness (posop)
  • Negative openness (negop)
rsaga.get.usage('ta_lighting', 5, env = env)
library path: C:\Users\b1066081\Desktop\SAGA-7~1.0_X\tools\
library name: ta_lighting
library     : ta_lighting
Usage: saga_cmd ta_lighting 5 [-DEM <str>] [-POS <str>] [-NEG <str>] [-RADIUS <double>] [-METHOD <str>] [-DLEVEL <double>] [-NDIRS <num>] [-UNIT <str>] [-NADIR <str>]
  -DEM:<str>        Elevation
    Grid (input)
  -POS:<str>        Positive Openness
    Grid (output)
  -NEG:<str>        Negative Openness
    Grid (output)
  -RADIUS:<double>  Radial Limit
    Floating point
    Minimum: 0.000000
    Default: 10000.000000
  -METHOD:<str>     Method
    Choice
    Available Choices:
    [0] multi scale
    [1] line tracing
    Default: 1
  -DLEVEL:<double>  Multi Scale Factor
    Floating point
    Minimum: 1.250000
    Default: 3.000000
  -NDIRS:<num>      Number of Sectors
    Integer
    Minimum: 2
    Default: 8
  -UNIT:<str>       Unit
    Choice
    Available Choices:
    [0] Radians
    [1] Degree
    Default: 0
  -NADIR:<str>      Difference from Nadir
    Boolean

And compute with defaults:

rsaga.geoprocessor(
  'ta_lighting', 5, 
  list(
    DEM='input_dem/dsm_sa1.sgrd', 
    POS='out_products/posop_sa1.sgrd',
    NEG='out_products/negop_sa1.sgrd',
    METHOD=0,
    DLEVEL=3
  ),
  env = env
)

Channels module

Vertical Distance to Channel Network

For this some previous steps are needed, see the complete code:

rsaga.get.usage('ta_hydrology', 0, env = env)
rsaga.geoprocessor(
  'ta_hydrology', 0, 
  list(
    ELEVATION='input_dem/dsm_sa1.sgrd', 
    FLOW='int_products/flowsa1.sgrd', 
    METHOD=0
  ),
  env = env
)

rsaga.get.usage('ta_channels', 0, env = env)
rsaga.geoprocessor(
  'ta_channels', 0, 
  list(
    ELEVATION='input_dem/dsm_sa1.sgrd', 
    CHNLNTWRK='int_products/chnet_sa1.sgrd', 
    INIT_GRID='int_products/flowsa1.sgrd', 
    INIT_VALUE=1000000
  ),
  env = env
)

rsaga.get.usage('ta_channels', 3, env = env)
rsaga.geoprocessor(
  'ta_channels', 3, 
  list(
    ELEVATION='input_dem/dsm_sa1.sgrd',
    CHANNELS='int_products/chnet_sa1.sgrd',
    DISTANCE='out_products/vdcnw_sa1.sgrd'
  ),
  env = env
)

Final steps

Useful function to convert to .tif extension.

files = list.files(path = "./out_products", pattern = "sa1*.*sdat", full.names = T)

convertToTiff <- function(x) {
    b = brick(x, crs = '+init=epsg:2193')
      b = projectRaster(b, crs = '+init=epsg:2193')
    extension(x) = '.tif'
    dt = dataType(b)
    writeRaster(b, filename = x, datatype = dt)
}

lapply(files, convertToTiff)

And fast, view the results:

files = list.files(path = "./out_products", pattern = "*sa1.*tif", full.names = T)
file_names = list.files(path = "./out_products", pattern = "*sa1.*tif")
library(stars)
library(purrr)
files_read = files %>% map(function(x) read_stars(x, proxy = T))

par(mfrow = c(5,6))
invisible(lapply(files_read, function(x) {
  plot(x, key.pos = NULL, reset = F, main = NULL)
}))

Names and SAGA code of generated output

list = read.csv('code_list.csv', sep = ';')
list

References

Eisank, C., Hölbling, D., & Friedl, B. (2014). How well do terrain objects derived from pre-event DEM spatially correspond to landslides? Geological Society of America Abstracts with Programs, 46(6), 105. Vancouver, Canada.